Leptonic widths of heavy quarkonia: S-Wave QCD/NRQCD matching coefficients for 
the electromagnetic vector annihilation current at 0(a s v 2 ) 
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We construct the S-wave part of the electromagnetic vector annihilation current to 0(a s v 2 ) on 
the lattice for heavy quarks whose dynamics are described by the NRQCD action, and v is the 
non-relativistic quark velocity inside the meson. The lattice vector current for QQ annihilation is 
expressed as a linear combination of lattice operators with quantum numbers L = 0, J p = 1~, and 
the coefficients are determined by matching this lattice current to the corresponding continuum 
current in QCD to 0(v 2 ) at one-loop. The annihilation channel gives a complex amplitude and 
a proper choice for the contours of integration is needed; a simple Wick rotation is not possible. 
In this way, and with a careful choice of subtraction functions in the numerical integration, the 
Coulomb-exchange and infrared singularities appearing in the amplitudes are successfully treated. 
The matching coefficients are given as a function of the heavy quark mass Ma in lattice units. An 
automated vertex generation program written in Python is employed, allowing us to use a realistic 
NRQCD action and an improved gluon lattice action. A change in the definition of either action 
is easily accommodated in this procedure. The final result, when combined with lattice simulation 
results, describes the electromagnetic decays of heavy quarkonia, notably the T meson. 

PACS numbers: 12.38.Bx, 12.38.Gc, 13.20.Gd 



I. INTRODUCTION 

Heavy quark states like the J/* [J H and T || 
mesons play a central role in the experimental study of 
the electroweak interactions. It is therefore important 
that we have reliable non-perturbative QCD predictions 
of their properties against which to compare. Lattice 
Monte Carlo simulations provide the only systematically 
improvable framework for such studies, but relativistic 
quark actions do not lend themselves very easily to lat- 
tice simulations of heavy quark dynamics; the Compton 
wavelengths of heavy quarks are small compared to cur- 
rently feasible lattice spacings. 

Fortunately the heavy quarks are much heavier than 
the hadronic scale A w 200 MeV, whilst their kinetic en- 
ergy is small (as demonstrated by the radial excitations 
of the mesons being much smaller than the ground state 
energy). This allows a non-relativistic description of the 
mesons using the NRQCD effective field theory 0, Q , us- 
ing the heavy quark velocity as the expansion parameter. 

Simply put, the goal of this paper is to provide match- 
ing coefficients that allow NRQCD matrix elements (cal- 
culated non-perturbatively in a lattice simulation) to be 
used to predict heavy quark phenomenology, in particu- 
lar the leptonic widths of the T mesons. 

More precisely, to obtain accurate results from a lat- 
tice simulation the QCD and NRQCD actions must be 
systematically improved to eliminate errors due to lat- 
tice artifacts, relativistic corrections and radiative effects. 
Both perturbative and non-perturbative methods exist to 
do this. A similar programme is needed for improvement 
of lattice operators and currents. In this paper we use 
perturbation theory to match matrix elements of the S- 



wave part of the vector QQ heavy-heavy electromagnetic 
annihilation current calculated on the lattice to the con- 
tinuum result, ensuring that the lattice results give the 
correct answer to 0{a s ) in the strong coupling constant 
and 0(v 2 ) in the velocity. This technique has already 
been used to improve the weak annihilation current for 
leptonic B-meson decay @, Q via the weak annihilation 
of a heavy quark Q and light anti-quark q. 

The QQ annihilation is more complicated than the 
weak Qq case. In the heavy-light case, we can exploit the 
crossing symmetry of the relativistic light quark action 
to match instead the weak heavy-light Qq form-factor. 
The amplitude for this is purely real and so the choice 
of integration contour for the temporal component of the 
momentum is straightforward (parallel to the imaginary 
axis). The NRQCD action lacks this crossing symmetry 
and so the time- like improved lattice vector current (rele- 
vant for annihilation) is not a priori related to its space- 
like counterpart (which determines the heavy quark form- 
factor). The amplitude for on-shell QQ annihilation is 
complex with a threshold for QQ scattering and has the 
additional complication that it contains a Coulomb sin- 
gularity. The calculation therefore entails a careful choice 
for the integration contours. For heavy quark velocities 
v > this does not correspond to the simple Wick ro- 
tation (generally with constant real part displacement) 
which suffices for the improvement of the form-factor. 

In addition, the Coulomb singularity gives rise to terms 
odd in v starting at C(w _1 ), and the integrand must be 
subtracted in a suitable way so that the numerical inte- 
gral along the contour that passes close to the singularity 
can be done accurately. None of these difficulties occur 
in the matching calculations for the space-like weak and 
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electromagnetic form factors involving heavy quarks. We 
choose to match the real part of a suitable linear combi- 
nation of the electromagnetic form- factors Fi(q 2 ), F2(q 2 ) 
for time-like q with q 2 = AM 2 (I + v 2 ). 

Earlier matchings of the vector annihilation current 
avoided these issue s by either being restricted to tree level 
@ or to v = [HI El- Neither is satisfactory: v 2 and a s 
are comparable at around 0.1 in the T system and failure 
to include both leads to strong discretisation errors in 
calculations of the leptonic width (T^ |. In addition, in 
flU only the simplest NRQCD action was used, keeping 
only terms to leading order in v. 

This study corrects this, using a gauge action that al- 
lows lattice matrix elements to be calculated using the 
state-of-the-art lattice QCD ensembles produced by the 
MILC collaboration. The NRQCD action is the same im- 
proved form used in recent studies, e.g. [H, [li|. When 
the matching coefficients calculated here are married to 
lattice NRQCD matrix elements, it will allow a determi- 
nation of the leptonic width that is correct to 0(10%) 
and of the ratio of the widths of the 2S and IS T states 
correct to within a few per cent. The size of these uncer- 
tainties matches those in the experimental measurements 
[T^ |. which justifies our one- loop, perturbative approach 
in the matching. 

The structure of the paper is as follows. In Sec. HT1 we 
describe the matching procedure. The continuum QCD 
matrix element analytic calculation is given in Sec. lIIII In 
Sec. lIVI we present the numerical calculation of the corre- 
sponding NRQCD matrix elements. The final matching 
coefficients are determined in Sec. [V] and discussed in 
Sec. lVIl In the appendices, we describe the tests we have 
applied in our calculation to ensure the correctness of the 
Feynman rules and of the loop integration, and also to 
establish the independence of the results on the gauge 
fixing and infrared regulator. 

A preliminary version of this work was presented in 
Ref. [l| . 



II. MATCHING S-WAVE DECAYS 

The leptonic width of a heavy quarkonium state of 
mass Mqq is given in terms of a QCD matrix element 
M M e by 

Tee = fUU? I MmE ! 2 e Q a om (!) 

QQ 

where eQ is the electric charge of the heavy quark and 
a cm the fine structure constant. The matrix element 
represents the probability of the heavy quarks meeting 
and annihilating, and in the simplest picture is repre- 
sented by a hydrogenic "wavefunction at the origin": 

|M ME | 2 * 1>H(0). 

To compare with the experimentally measured widths, 
we want to calculate this matrix element in continuum 
QCD, in a way that embodies all the non-perturbative 



dynamics. As explained in the introduction, we cannot 
do this directly and must instead use lattice NRQCD 
simulations. 

The problem is that we do not know a priori which 
NRQCD current we should use. Instead, we should con- 
sider a set of suitable currents and separately calculate 
the matrix elements of each one using the Monte Carlo 
generated ensemble. The true QCD matrix element is a 
linear combination of these, and this paper provides the 
necessary coefficients. 

We choose our NRQCD currents to be 

* = (2) 

where bold face symbols denote spatial 3- vectors and M 
is the heavy quark mass. 

To convert our non-perturbative lattice NRQCD cur- 
rent matrix elements into the corresponding QCD value, 
Mme = (0 | jQ CD | QQ), we need matching coefficients a* 
such that 

(o\j QCD \QQ) =£a i (0|J i |QQ) . (3) 

i 

In this paper we fix them. 

When we calculate the matrix elements of Ji in the 
simulation, the mass M will of course be replaced by 
a number. We may choose it to be the bare mass or 
(less usually) the renormalised value. The a t will differ 
accordingly, so we will give separate results for the bare 
and renormalised cases. 

Our matching method is summarised as follows: the 
NRQCD matrix elements each depend differently on 
the heavy quark velocity (at tree level, for instance, 
(0\Ji \ QQ) oc v 21 ). By choosing the a; appropriately we 
can match the QCD velocity dependence order by order 
in v 2 . We make this choice perturbatively, performing 
the velocity matching at each order in a s in turn 15]. 

We start by expanding the currents and matching co- 
efficients as power series in a s : 

En (n) 

n 

(o\j\qq) = Y, a s(°\ J \QQ) in) ■ ( 4 ) 

n 

The superscript (n) denotes the 0(a") perturbative con- 
tribution. 

Working in the Breit frame (where the decaying meson 
is stationary), we take the Euclidean four momentum of 
the quark and antiquark to be 

jf = (iE Q ,±p) , p= (0,0, aMv) . (5) 

The dimensionless expansion parameter is v. Although 
we refer to it as the heavy quark velocity, it is related to 
the true velocity by v = /3/y/l — {3 2 . 
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We treat the heavy quarks as being on shell. This is 
exact, even though we might expect off shell contribu- 
tions at 0(a 2 ). By using the equations of motion, the 
contributions from these within a bound state are seen to 
vanish at all but a subset of spacctime points of measure 
zero EEIT3I. 



A. Matching at tree level 

The matching at tree level is essentially trivial. Using 
the standard Dirac representation for the 7-matrices in 
terms of Pauli er-matrices: 



7° = 











u 


X) 




-a* 



(6) 




(7) 



where tp and x are the standard Pauli spinors for quarks 
and antiquarks, respectively. We have chosen the non- 
relativistic normalisation for consistency with NRQCD, 
since the Foldy-Wouthuysen-Tani transformation [3, 0, 
[H SJ is unitary. 

In terms of these Pauli spinors, the relevant Dirac 
tensor components of the non-relativistic expansion of 

the tree-level matrix element (0 I jQ CD ;^| , u 1 
v(—p)^f^u(p) are: 



v(-ph°u(p) 

v{-p)iu(p) 



= 



i #'2 M\ 



v(-p)- 



ia i0 E 



M 



= fi{v 2 ) xW , 

MP) = + 

= / 2 (« 2 )xVV. 



(8) 



where we have averaged over spatial directions for S-wave 
decays [§]. 

The trcc-lcvcl matching coefhcients must satisfy the 
leading order term in Eq. ([3]): 



(0|J QCD |QQ> (O) = 5>f (0\Ji\QQ) 



AO) 



(9) 



The expansions in powers of v are 
fi(v 



2 ' = l-\v 2 + \v* + 0{J), 



h{v 2 ) = l + \v 2 ~^v^O{v"). (10) 
Using Eq. (0), the tree level NRQCD matrix elements 
can be written as 



<0|J,|0Q) (0) =9i(v) xW 
The tree level velocity dependence is 
g (v) = 1 



(11) 



92{v) 



4 o f aMv\ 
-o {—) 

2 ( aMv\ 



(aM) 2 
4 



(aM) 4 



4 sin 



^— — J - sin (aMv) 



(12) 



such that gi(v) — (—v 2 ) 1 at lowest order in v. 

A term by term comparison of these expansions with 
that of /1 yields 



,(0) 



,(0) 



(aMY 
72 



(13) 



B. Matching at one- loop order 

To match at one-loop order, we need to calculate the 
one-loop QCD and NRQCD corrections to the quark- 
antiquark annihilation vertex. The QCD corrections con- 
sist of both self-energy insertions on the external legs and 
a vertex correction, and for the case of a quark-antiquark 
vertex can be written as 



QCD 



FIG. 1: One-loop corrections to the quark-antiquark annihi- 
lation current in QCD 






(0|J QCD |QQ) (1) = F^(4E 2 )v(-p) 7 u(p)+iF^(4E 2 ) v{-p)qu{p) 
F^iAE^Mv 2 ) + F^(4E 2 )f 2 (v 2 )} X W 



(14) 
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where q\ = <7i U q v /M and F^ are the 0(a s ) contributions to the vertex structure functions. We note that while after 

renormalisation F$p is UV-finite because of the Ward identity, it contains IR divergences. These infrared divergences, 
however, are the same as those that arise in NRQCD, since the low-energy behaviour of the two theories is the same. 
The 0(a s ) matching condition from Eq. ((3]) is then 



(o\Ji\QQ) 



(0) 



F«(4£: 2 )A(^) + F^(4E 2 )f 2 (v 2 )] X W " E a P <° ^ 

i 

(Iqcd - ^nrqcd) X^°"0 



(i) 



(15) 



The infrared divergences cancel between the first two 
terms, leaving an IR- and UV-finite expression that can 
be evaluated numerically. We opt to project out the a 2 
component and use the tree-level expectation values of 
the NRQCD operators Ji as our basis functions to fit 
the difference between the QCD and NRQCD one-loop 
results and determine 



E fl i = J Q CD ~~ /] 



NRQCD 



(16) 



To match to 0(v 2 ) in this calculation i runs from to 1 
only. 



III. CONTINUUM QCD CALCULATION 

To evaluate Iqcd analytically, we must regulate the in- 
frared Coulomb divergence in the Feynman integrals. To 
avoid the complications of twisted boundary conditions, 
we introduce a gluon mass fi and use the gauge invariant 
Stiickelberg propagator for the massive vector field (see 
Sec. (3-2-3) of Ref. H|): 



9fiv k^kjj / /i 2 
k 2 — u? + ie 



k^k^/ /x 2 



k 2 -/i 2 /A + 



ie 



(17) 



where A is the gauge fixing parameter. 

The one-loop QCD contribution is given by the sum 
of the Feynman diagrams shown in Fig. [TJ The two left- 
most rescale the tree-level element by the quark wave 
function renormalisation constant Z . The rightmost dia- 
gram is the one-loop vertex correction. The full one-loop 
vertex function is a rather formidable-looking expression 
(HJ . We know from the Ward identity, however, that the 
vertex function must take the form 



u(p')r M u(p) = u(p') 



—F 2 (q 2 )a^ 



u{p) 



when sandwiched between on-shell spinors, where q = 
p — p' is the gluon momentum flowing out of the vertex. 
We also know from the Ward identity that Z^ 1 = -Fi(O), 
so that we can renormalise the vertex function order by 
order by subtracting from F\(q 2 ) its value at zero gluon 
momentum to obtain the rcnormalised structure function 



F[ n \q 2 ) 



F[ n) {$) . (19) 



This amounts to including the effects of the first two 
diagrams, with which we will therefore no longer concern 
ourselves. 

For the case of quark- antiquark annihilation, we then 
have 



<0|J m qcd |QQ) (1) =«(-p) 



F; 



iE 



i ^E 2 ) ltl + —F^(AE 2 )a^ 



u(p) 



(20) 



or, in terms of Pauli spinors 



J^ CD QQ 



(i) 



X W \F[^ R {AE 2 )h{v 2 ) + F^{AE 2 )f 2 (v 2 ) 



(21) 



To compute F\ and F2 without resorting to the Feyn- 
man or Schwinger parameter representations (which are 
not available for NRQCD because the denominators are 
not quadratic), we employ a number of techniques. A 
discussion of these will be useful later. 



Since the decomposition of the vertex function into 
form factors stated above is only valid between on-shell 
spinors, we put it between the appropriate on-shell pro- 
jectors 
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(f + M)T^p\p){i, + M) = {{,' + M) 



{i> + M) 



(22) 



r 



where the appropriate on-shell momenta for an incoming 
quark-antiquark pair are given by Pfl = (E, p) and p' p — 

{-E,p) with E = ^/M 2 +p 2 . 

Contracting the above equation with either {p + p'Y 
or 7^, and taking the trace of both sides, we obtain two 
equations for F\ and F 2 : 

A = {P + P ' Y Tr((/ + M)Y ll { P ',p){i > + M)) 



2M 
(p + p'f 



(4M 2 F 1 (q 2 )-q 2 F 2 (q 2 )) 



2M 2 

B = Tr( 7 ^' + M)iyp',p)(^ + M)) 
= 4(2M 2 + q 2 )F 1 (q 2 )-6q 2 F 2 (q 2 ) , 



(23) 



with solutions 
Fl{q 2 ) = 
F 2 (q 2 ) = 



4(4M 2 - q 2 ) 
2M 2 (2M 2 + q 
q 2 (p + p') 2 
A 



12M 2 
(P + P') 2 



A-B 



B 



{p + p'f 2(2M 2 + q 2 ) 



(24) 



At the one-loop level, the vertex function is given by the 
integral expression (in Feynman gauge) 



AttCo 



d 4 k 7p (/' + M) 7m (/ + M) 7 p 

(2tt) 4 (k 2 - fi 2 )(l' 2 - M 2 )(l 2 - M 2 ) 

(25) 

where we have defined the loop momenta I = k + p and 
I' = k+p' and introduced a gluon mass fi as an infrared 
regulator. After performing the manipulations outlined 
above, this vertex function leads to 



) 4 (p + p') 2 {k 2 - M 2 )(Z' 2 - M 2 ){1 2 - M 2 ) 



6M 



F^(q 2 ) = 4ttC 2 



2M 2 l-l' + 2M 2 (p + p')-(l + l') f (/)2 

+2p • Z'p' -l-M 4 - M 2 p ■ p'] , 

d A k 2M 2 {2M 2 + q 2 ) 



( p + p').l(p + p')-l f 



{2ir) A q 2 (p + p'f{k 2 - M 2 )(Z 2 - M 2 )(l' 2 - M 2 ) 

2 



4 ((p + p') ■ (I + I') + I ■ I' — M 2 ) 



(P + P') 



(p + p>).l(p + p>).l> 



2M 2 + q 2 



(M 2 l ■ I' - 2p ■ I' p' ■ I + M 2 (p + p') ■ (l + l')+p-p' - 2M 2 ) 



(26) 



for the structure functions. 



In the physical limit /j, — > 0, the one- loop structure 
functions are of course well-known analytically, since they 
are just the QED structure functions multiplied by the 



group-theoretic factor C 2 = 4/3: 



f h % 2 ) = 



F^\q 2 ) 



log-^ + l) (0 cot 0-1) 

/ 

Jo 



g 2 C 2 

4tt 2 

rO/2 

+2 cot 9 I otan 
g 2 C 2 6 



9 
— tan — 
4 2 



8tt 2 sin0 



(27) 



6 



where 



6 = 2arcsin(£;/M) 



(28) 



We have compared our numerical evaluation of the struc- 
ture functions in both the form factor and annihilation 
channels with their analytical values and have found ex- 
cellent agreement. Especially, we were able to replicate 
the infrared divergence by varying our gluon mass \x. Re- 
solving the 1/v Coulomb singularity in the annihilation 
channel requires special care. To avoid a contamination 
of the low-v behaviour by the gluon mass fi, which acts 
as a cut-off on the v dependence by limiting the momen- 
tum of the exchanged gluon, we have scaled fi with v, 
and then were able to observe the correct Coulomb sin- 
gularity behaviour in the infrared finite part of F± . 



A. Wick rotation 

In doing these calculations, we must be careful how 
we Wick rotate our integration contour. In the quark- 
antiquark annihilation channel, the poles of the inte- 
grands in the complex fco plane are located as shown in 
Fig. [2] For (fe + p) 2 > p 2 , the poles are all located 
second and fourth quadrants of the Argand diagram for 
ko, and the usual Wick rotation of the integration con- 
tour is possible as in Fig. [2^. When (fe + p) 2 < p 2 , 
the fermionic poles cross the imaginary ko axis and we 
need to be more careful and deform the contour as per 
Fig. This choice of contour is, however, impractical. 
The short piece of the contour running along the real 
axis is by far the most dominant contribution to the in- 
tegral. We will estimate the value of the integral using 
Monte Carlo methods. To get this contribution correctly, 
we need to sample all three-momenta along the contour 
with comparable weights. We therefore use the equiva- 
lent contour shown in Fig. [2t, which works much more 
efficiently. 

In this triple contour case, we choose the outlying con- 
tours to be midway between the gluonic and fermionic 
poles. The Stiickelberg gluon propagator has two poles: 
one associated with the physical gluon mass, and a sec- 
ond at // 2 /A. To avoid possible numerical instabilities, 
we use the smaller of the two to fix the position of the 
outer two contours. 

Note that if we work with v = as in Ref. we 
can always Wick rotate as per Fig. [2ji. For non-zero v, 
however, it is important to note that the choice of an ap- 
propriate contour is essential to obtain the correct result: 
with a naive standard Wick-rotation, the structure func- 
tions obtained would not even be Lorentz-invariant. We 
have explicitly checked that our choice of contours leads 
to structure functions that are invariant under a Lorentz 
boost. Another point to note is that even in the quark 
form factor channel at spacelike q 2 , where the quark poles 
do not cross each other, a standard Wick rotation about 
the origin is not correct, and the rotated contour has to 
be shifted along the real axis by an amount depending on 



the kinematic frame, in order to pass between the poles 
and pick up the correct result. 



IV. LATTICE NRQCD CALCULATION 

In this section we describe the perturbative calculation 
using the lattice NRQCD action. 



1. The NRQCD Action 

The NRQCD action we consider is the same as Gulez 
et al. [HI], and also the same as has been used in recent 
simulations [HI, [2t| (although there is a typographical 
error in the description in the latter [24|): 



Snrqcd = ^2 ~ ^ ( 1 

x,t ^ 

aH 



xf/](l 



2?; 



iSH 

1 - 



1 - 



aHo 
2n 



aSH\ , s 

—J V> . (29) 



The i$ field is understood to be be located at (t, x) , with 
the position of ip on the timeslice t — 1 fixed by gauge 
invariance. Other than consistency with previous work, 
there are no strong arguments for the relative ordering 
of the kinetic and interaction terms in the action. The 
ordering here differs from, for instance, Ref. 

The leading kinetic term is 



H a = — 



A 2 
2aM 



(30) 



where M is the bare mass and n is a stability parame- 
ter for the nonrelativistic evolution equation, that must 
fulfil the condition n > 3/ (2aM) for the time reversal 
symmetric evolution equation [rll25l]. Gluonic corrections 
decrease the lower bound on n to just above 1 / (aM) [26| . 

The time-reversal symmetric splitting of the Hq oper- 
ator either side of the temporal link Q is designed to 
mimic the full time evolution due to Hq along a tempo- 
ral lattice spacing in a way that avoids the well-known 
instability in the discretisation of parabolic differential 
equations (see, for instance, Sec. 19.2 of Ref. [13])- In 
this way, the time step in the evolution equation is small 
enough to allow the highest momentum modes in the the- 
ory to come into equilibrium, whilst avoiding the need for 
a very small lattice spacing which makes the theory too 
expensive to simulate. 

The interaction term corrects for relativistic and dis- 
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FIG. 2: Locations of the poles and choice of integration contour for the (Minkowski metric) ko integration in QCD quark- 
antiquark annihilation. The solid circles represent poles in the gluon propagator, and the open circles the fermionic poles for 
various Ifcl. 



crctisation effects: 



a6H 



-ci 



(A 



(2h2 



C 2 



V -E-E V 



8(aM) 6 8(aNiy 
C3 1 , <2 <r • (V x E - E x V) 



-Ci 



+c 5 



8(aM) 
1 

2(oM) 
AW 



B 



C6 



(A 



(2)\2 



24(aM) 16n(aM) 



(31) 



We note that improved derivatives are used in the term 
proportional to C3 and that the improved field strength 




vertex 



T X) T 

^earlobe \ x2 1 bubble 




1 




tadpole 



FIG. 3: One-loop corrections to the self energy and annihila- 
tion current in NRQCD. The gluons in these diagrams can be 
temporal as well as spatial. The solid (blue) circles represent 
the current in Eq. <(3j . The open (red) circle represents the 
contribution from tadpole improvement of the current. "x2" 
denotes a similar diagram on the outgoing fermion line. 



has not been rendered explicitly traceless. The terms 
proportional to Cj for i = 1...4 provide relativistic correc- 
tions to 0(Mv 4 ) [1,[28| and represent the relativistic cor- 
rection to the kinetic energy, the nonabelian analogue of 
the Darwin term, the spin-dependent interactions lead- 
ing to spin-orbit couplings and the quark chromomag- 
netic moment, respectively. The final two terms remove 
the leading order discretisation error. All terms are un- 
derstood to be tadpole improved. We use the tree level 
values Ci = 1, as in Refs. [H, 03 ■ Other than tadpole im- 
provement, we do not consider the effects of radiatively 
correcting the Cj. 

We obtain the Feynman rules for the NRQCD and 
gauge actions using an automated procedure [301 ] , as out- 
lined in Appendix [AJ We also detail there the tests we 
have carried out to ensure that the Feynman rule expres- 
sions are correct, and the techniques we employ to speed 
up their evaluation for specific momenta. 



2. The lattice gauge action 

To maintain compatibility with the MILC collabora- 
tion simulations, we use the Symanzik improved gauge 
action 



+0(a s ) 



(32) 



where P, R are lxl and 2x1 Wilson loops respectively. 
0{a s ) denotes possible radiative and tadpole improve- 
ment of the action. As discussed later, these terms will 
not contribute at one-loop. 

The inverse lattice Stuckelberg propagator is 



r^(fc) = V^(k) + (a^yS^ + \k^k v 



(33) 



where the two-point function V^" depends on the action 
chosen. Gauge invariance requires the gauge fixing term 
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TABLE I: Diagrams contributing to matching calculation, 
precision of the number. 



Where no statistical error is given, it is smaller than the quoted 



aM 


V 


(IqCU~ (-^ vertex — -fin) ^out 


h 


^earlobe 


-^bubble 


^tadpole 


^QCD — ^NRQCD 






/odd) 













4.0 













0.02155 (1) 


-0.06695 (3) 


0.04688 






0.03 


-0.8511 


-1.3258 (28) 


1.2176 (5) 


1.8226 (9) 


0.0216 


-0.0668 


0.0468 


-0.1315 (30) 




0.07 


-0.8611 


-1.3112 (39) 


1.2262 (5) 


1.8212 (9) 


0.0220 


-0.0661 


0.0463 


-0.1452 (40) 




0.10 


-0.8738 


-1.2979 (44) 


1.2370 (5) 


1.8247 (9) 


0.0224 


-0.0652 


0.0456 


-0.1615 (45) 




0.15 


-0.9044 


-1.2636 (55) 


1.2623 (5) 


1.8288 (9) 


0.0234 


-0.0631 


0.0441 


-0.2004 (56) 


2.8 













0.05203 (2) 


-0.13678 (5) 


0.09566 






0.03 


-0.8511 


-1.3689 (21) 


0.9240 (4) 


1.6225 (9) 


0.0521 


-0.1365 


0.0956 


-0.1736 (23) 




0.07 


-0.8611 


-1.3678 (29) 


0.9468 (4) 


1.6241 (9) 


0.0526 


-0.1359 


0.0951 


-0.1809 (31) 




0.10 


-0.8738 


-1.3659 (33) 


0.9550 (4) 


1.6249 (8) 


0.0532 


-0.1349 


0.0944 


-0.1862 (35) 




0.15 


-0.9044 


-1.3622 (39) 


0.9726 (4) 


1.6266 (9) 


0.0545 


-0.1327 


0.0929 


-0.2007 (40) 


1.95 













0.12394 (3) 


-0.28210 (6) 


0.19724 






0.03 


-0.8511 


-1.3636 (15) 


0.7415 (2) 


1.3505 (8) 


0.1241 


-0.2820 (1) 


0.1971 


-0.1356 (17) 




0.07 


-0.8611 


-1.3669 (23) 


0.7459 (2) 


1.3486 (9) 


0.1246 


-0.2812 (1) 


0.1966 


-0.1357 (24) 




0.10 


-0.8738 


-1.3734 (26) 


0.7503 (2) 


1.3503 (9) 


0.1254 


-0.2804 (1) 


0.1960 


-0.1379 (27) 




0.15 


-0.9044 


-1.3979 (84) 


0.7634 (4) 


1.3537 (14) 


0.1271 (1) 


-0.2780 (1) 


0.1944 


-0.1319 (86) 


1.0 













0.5093 (2) 


-1.0720 (4) 


0.75000 






0.07 


-0.8611 


-1.3681 (13) 


0.5004 (2) 


0.4335 (10) 


0.5103 (2) 


-1.0713 (4) 


0.7494 


0.4040 (17) 




0.10 


-0.8738 


-1.3956 (16) 


0.5034 (2) 


0.4323 (10) 


0.5110 (2) 


-1.0701 (4) 


0.7488 


0.4044 (19) 




0.15 


-0.9044 


-1.4146 (18) 


0.5127 (2) 


0.4323 (10) 


0.5132 (2) 


-1.0676 (4) 


0.7472 


0.4005 (21) 



to be constructed from lattice momentum vector = 
2sin(afc M /2), so the Feynman gauge (A = 1) propagator 
is only diagonal for the Wilson gauge action (for which 
the gluon two-point function is = k p k p 5 111 ' — k^k u ). 
Note that A = 1/a in the notation of Ref. [3l|. As we 
do not consider Landau gauge here, the inverse propaga- 
tor is directly invertible and we do not need to use an 
intermediate gauge. 



3. Annihilation currents and radiative improvements 

We use lattice NRQCD annihilation currents that are 
the naive discretisations of Eq. Q : 



J = ^ZXI^OL 



^ Ax (aM) 2 

x;i=l y ' 

x (Ui{x)i) x+i + U}(x - i)i) x - % - 2^) (34) 

and the links in J\ are understood to be tadpole im- 
proved. Removing the mean field, "tadpole" contribu- 
tions improves the convergence of lattice perturbation 
theory markedly (32J. Operationally, this is done by di- 
viding every gauge link U in the action by a factor of uq. 
Common definitions for uq are that it is the mean link in 
Landau gauge or the fourth root of the mean plaquette. 
We use the former, expanding the link perturbatively as 



(2) 



it = l-a s % 
as used Ref. 11 



with u. 



(2) 




0.750 from Ref. \m and 



Tadpole improvement of the NRQCD action does not 
contribute to our calculation, as the fermion wavefunc- 
tion rcnormalisation has no tadpole correction for the 
time reversal symmetric form of the NRQCD action (the 
argument mirrors the mean field analysis in Ref. @ ) . As 
discussed before, we do not consider any further radiative 
improvements of the NRQCD action. 

Tadpole and other radiative improvements of the gauge 
action also do not contribute to the matching calculation. 
The leading order effect of these is an 0(a s ) insertion in 
the gluon propagator. As there are no external gluons 
in our calculation, such insertions will only contribute at 
two loops and above. 

The only effect of tadpole improvement comes from the 
current J\ , and its contribution to /nrqcd can be easily 
calculated: 



'tadpole 



lufaf 
(aM) 2 

(2) (0) 



cos pi 



(aM)< 



+ v 2 + 0(v*)) .(35) 



The only other possible source of radiative corrections 
comes from the mass used in Eq. ^ when we calcu- 
late the non-perturbative NRQCD matrix elements in 
the Monte Carlo lattice simulation. If the number M 
used in the simulation is the renormalised heavy quark 
mass, there is no further correction. If the number for 
the bare mass is used instead, the renormalised mass 
is ZmM and we should divide the matching coefficient 
ttj by (Z M ) 2 ' 1 . In this study, that amounts to shifting 



9„(0) 7 (1) 



We calculate the multiplica- 



tive mass renormalisation factor in Appendix [Bl and will 
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present our results for the matching coefficients both with 
and without the shift. 



A. Calculating the vertex corrections 

The one-loop diagrams contributing to the quark- 
antiquark annihilation amplitude in NRQCD are shown 
in Fig. [3l The NRQCD fco integrals are around the unit 
circle in the e zfe ° a -plane. The quark poles sometimes cross 
the unit circle (just as the they crossed the imaginary axis 
in the QCD integrals), so we scale the circle of integration 
to avoid them and adopt a similar triple-contour strat- 
egy: integrating along three appropriately scaled concen- 
tric circles when the poles cross each other, and along the 
unit circle otherwise. 

The gi(v) are all even functions of v, but /qcd and 
^nrqcd both contain odd powers. We must assure our- 
selves that these exactly cancel in Eq. (|T5|) . The argu- 
ment is that NRQCD is an effective theory of QCD which 
can be systematically improved to reproduce all features 
of QCD, including the odd powers. There are, however, 
no S-wave operators containing odd powers of v that we 
could use in the improvement. The odd powers must 
therefore cancel exactly in Eq. (fT"5|) . This is not entirely 
surprising given that the odd powers arise from an even 
polynomial in v multiplied by the 1 /v Coulomb 1R diver- 
gence, and we know that NRQCD must reproduce the IR 
physics exactly. Nonetheless, it is worth examining the 
cancellation in more detail. 

Consider the power-expansion of the QCD expression 

F W- R {AE 2 )h{v 2 ) + F^{AE 2 )h{v 2 ) (36) 

with Fip' R defined in Eq. (fT9]) . Using the analytic results 
given above, we see that both fi(v 2 ) and f2(v 2 ) contain 
only even powers of v. Any odd powers in the expansion 
must therefore come fromF 1 (1) " R (4 J B 2 ) or F 2 (1) (4E 2 ). The 
analytical evaluation of the structure functions shows us 
that (4E 2 ) contains odd powers in v only in its imag- 
inary part, so the odd powers in the final answer must 
come from F[ X) ' R {AE 2 ). 

On the NRQCD side, we know that any odd powers 
in v must come from the quark pole giving rise to the 
Coulomb singularity, since the residues in the /co-plane 
of all other poles can be expanded in powers of v 2 . The 
odd powers therefore originate exclusively from integrals 
of the form 

r d*k (k 2 r( P 2 r ( r , 

J (2tt) 3 fc 2 (fc 2 + 2k-p + ie) ' 1 ' 

A careful analysis of these shows that only those integrals 
with a — contribute to the real part, whereas the others 
(which arc UV-divergent in the continuum) contribute 
only to the imaginary part. Since 

f d 3 k 1 1 

J (2tt) 3 fc 2 (fc 2 + 2k-p + ie) 16\p\ 1 ' 



the only odd powers of v in the NRQCD result will come 
from multiplying powers of p 2 in the numerator with the 
Coulomb singularity. Expanding these, we find the same 
coefficients multiplying each odd power of v as in the 
above QCD result. 

We note that to obtain correct results to order v 2n we 
have to use the correctly matched 0(v 2n ) tree- level anni- 
hilation operator. We must also use 0(v 2n ) quark-gluon 
vertices in the diagram involving spatial gluons and the 
0(v 2n+2 ) quark propagator (the expansion of the lat- 
ter around the Coulomb singularity pole gives an 0(v 2n ) 
contribution) . 

In summary, then, matching at tree- level to 0(v 2p ) 
guarantees the cancellation of the odd powers at one-loop 
level to C(w 2 p- 1 ). 

We will estimate the NRQCD loop integrals stochas- 
tically u sing the adaptive Monte Carlo package called 
VEGAS [2II2I (see Sec. 7.8 of Ref. ^ for further dis- 
cussion). These estimates of the NRQCD integrals will 
only converge if the integrands are both finite and rela- 
tively smooth. Both Iqcb an d /nrqcd have an infrared 
Coulomb divergence. Although these are formally regu- 
lated by the gluon mass, the integrands are still sharply 
peaked, leading to unacceptably slow convergence of the 
numerical integration. 

As we have discussed, all odd powers of v cancel point- 
wise in the difference of the two integrands, -Zqcd — 
-^nrqcd i leaving a smooth integrand. The obvious strat- 
egy is to numerically estimate the difference as a single 
integral, remembering that the NRQCD integrand is only 
defined inside the finite Brillouin zone. 

Direct subtraction has problems. The NRQCD inte- 
grand is quite complicated and time-consuming to eval- 
uate for given momenta. This limits the number of inte- 
gration points that VEGAS can consider in a set time. 
Conversely, the QCD integrand needs a large number of 
points to accurately estimate the integral: the terms like 
l/(p + p') 2 in Eq. (|2"6")l give rise both to an apparently 
UV divergent contribution to the 1 /v Coulomb singular- 
ity and l/v 2 term in the result. These terms, however, 
come with a factor of cos 9 from the scalar products with 
(p + p'), and thus vanish only after integration over all 
spatial angles. 

If we directly subtract the integrands, we arrive at a 
function that is both expensive to evaluate and needs 
many integration points to converge. To get round this, 
we use an analytic form of the QCD structure functions 
and only evaluate the NRQCD integrals numerically. For 
the latter we need to smooth out the regulated l/v in- 
frared divergence by subtracting an integrand with the 
same low-momentum structure. Fortunately, we can still 
cancel all the odd powers of v from the NRQCD inte- 
grand by multiplying the integral to be subtracted by an 
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appropriate function of v 2 : 

Ah(v 2 ) f d A k 



7 odd = Im { - 



x ik 



x [ ik + 



k 2 + 2kp 

2M 

k 2 + 2k ■ p 

2M 



h(v 2 ) 
12v 



(39) 



with 



h(v 2 ) 



[l + 2v 2 ) (l + 2VTT^) 
3(1 + v 2 ) 



This is certainly sufficient for the low powers of v 2 in 
which we are interested here. By comparing the respec- 
tive power-series expansions term by term, it can easily 
be seen that the odd powers of v are the same as in the 
QCD result. 

To evaluate Eq. (|16[) we therefore take the difference 
of (/qcd - /odd) calculated analytically and (Jnrqcd - 
/odd) estimated numerically. Both expressions are even 
power series in v, and the subtracted NRQCD integrand 
is now sufficiently smooth that no change of variable in 
the momentum coordinate, designed to "squash" many 
evaluation points onto the contour in the neighbourhood 
of the pole [36j |. is required. It is convenient to split 
J dd into two integration regions, within (/ n ) and outside 
(7 out ) the NRQCD Brillouin zone |fc^| < ?r/(aM). 

^vertex and Iz separately have infrared "cutting" di- 
vergences that cancel in their sum. Although the diver- 
gences are regulated by the gluon mass, by evaluating 
/vertex and Iz together we would have a smoother inte- 
grand for VEGAS. We meet the same problem as before, 
however: /vertex has a relatively cheap integrand but the 
VEGAS estimates are slow to converge. Iz converges 
quickly, but taking derivatives of Feynman rules makes 
the integrand expensive to evaluate. Therefore, we cal- 
culate the NRQCD integrals in Fig. [3] separately using 
VEGAS, choosing the number of integration points to 
give comparable statistical accuracy in the results. 

The final calculation is then made up of 



/QCD — /NRQCD — (IqCD ~ /odd 

= (/QCD - /odd) — 

([/vertex /n] /out 
+/bubblc + /tadpole) • 



RESULTS 



(^NRQCD — /odd) 
h 



I 



carlobc 

(41) 



In this paper we present results for four choices of 
heavy quark mass: aM = 4.0, 2.8, 1.95 and 1.0. The 



TABLE II: The matching coefficients, as a function of the 
renormalised heavy quark mass, for the leptonic width (a,i) 

and leptonic width ratio {hi). Note that do = 1, o? = of 1 = ^, 
and that there is no subtraction to prevent mixing down. 



Ma 



.CO 



,(()) 



4.0 2 -0.1288 (27) -3.29 (29) 

2.8 2 -0.1732 (21) -1.27 (21) 

1.95 2 -0.1358 (16) -0.02 (16) 

1.0 4 0.4056 (20) -0.22 (16) 



-3.27 (30) -0.09722 

-1.24 (22) 0.01611 

0.00 (17) 0.07219 

-0.29 (17) 0.11111 



TABLE III: The matching coefficients, as a function of the 
heavy quark mass, for the leptonic width (en) and lep- 



tonic width ratio (bi). Note that do 



a? 



(0) 



and 



(40) that there is no subtraction to prevent mixing down. 



Ma 



1 o 



,,(!> 



4.0 2 -0.1288 (27) -3.32 (29) 

2.8 2 -0.1732 (21) -1.35 (22) 

1.95 2 -0.1358 (16) -0.16 (16) 

1.0 4 0.4056 (20) -0.50 (16) 



-3.30 (30) -0.09722 

-1.32 (22) 0.01611 

-0.14 (17) 0.07219 

-0.56 (17) 0.11111 



first three represent the 6-quark mass on the MILC im- 
proved staggered ensembles with a ~ 0.09 fm ( "fine" ) , 
0.12 fm ("coarse") and 0.17 fm ("super-coarse") 
Mass aM = 1 represents the charm quark mass on the 
super-coarse lattices. In agreement with Ref. [lH, we use 
n = 2 for all masses except aM = 1.0, where n = 4. 

We choose IR gluon mass {an) 2 = 10~ 4 and use Feyn- 
man gauge A = 1. In Appendix [C] we show that our 
results do not depend on either of these choices. We also 
compare with relevant existing results in the literature 
for v = 0. 

The NRQCD diagrams were evaluated for a range of 
velocities from v = 0.03 to v = 0.15. In addition, we 
evaluated /bubble and / ea riobe at v = 0. The results are 
shown in Table |TJ We extracted the matching parameters 
using a linear fit as per Eq. (JTHJ). The matching coeffi- 
cients a^ and a^ are given in Tables ITT1 and 1 1 1 1 1 Results 
from the former are to be used when the number for the 
renormalised heavy quark mass is used to construct the 
currents in Eq. (f3"4")l . Results from the latter are to be 
used if the bare mass is instead employed. 

The results in the renormalised mass case are shown 
graphically in Fig. 0J We note that for smaller masses, 

the coefficients of the J\ current, a± and 6j , are much 
smaller when renormalised masses are used. 

We have checked that the fits are not biased by higher 
terms in the velocity expansion. Note that whereas —Iz, 

/vertex 

and —/bubble reduce monotonically as aM is in- 
creased, /out, -/carlobc and -/tadpole grow. Given that 
the result of combining these will depend on (aM) 2 , 
(aM) 4 and l/(aM) 2 , it is not surprising that the match- 
ing coefficients do not vary monotonically with the heavy 
quark mass. 

Our computations of these diagrams have been per- 
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formed on the SunFire Galaxy-class supercomputer at 
the Cambridge-Cranficld High Performance Computing 
Facility using an implementation of the VEGAS algo- 
rithm adapted to parallel computers using MPI (Message 
Passing Interface). 



A. Mixing downwards 

Whilst at tree level matrix elements of J\ contribute 
only at 0(v 2n ), at higher loop orders there will be con- 
tributions at lower orders of v 2 . We call this "mixing 
down". In the case of Ji, the integrals I e ariobe, /bubble 
and /tadpole are only weakly momentum dependent and 
Jvertex also makes a contribution at v = 0. 

This is theoretically inconvenient as we must redo all 
previous calculations when we improve the current to 
higher orders of v 2 and cannot easily compare the new 
numbers with the old to check for consistency. We can 
get around this by introducing subtracted currents to 
prevent this downward mixing of currents. Although not 
essential for lattice Monte Carlo calculations, subtracted 
currents are also useful here as they make the conver- 
gence of the double series in a s and v 2 in Eqn. ([3]) most 
explicit. Thus we can expect that the matrix element of 
the subtracted J\ will vary as v 2 (to some order in a 

We define the subtracted currents as J = z. 
the coefficients z.y are chosen to prevent this downward 
mixing of currents at all radiative orders: 



zjJj, where 



(0\Ji\QQ) 



,.2i 



(42) 



At tree level z^ = 
z\a = for j > 



Sij. At higher loop level we set 
i, as we are only concerned with 



preventing downward mixing. For 0(v 2 ) matching, the 
only non-trivial element is z^ , fixed by 

z[f <0|J,|QQ) (1) +zW<0|J 3 |QQ) (O) 
- ^--(OlJilQQ)' 1 ' 



= 0, 



-10 



v=0 



(43) 



Note that in this calculation we consider only 7bu 



bble, 

^cariobc ^tadpole and I V crtcx, all at v = 0. Data for z[q are 
given in Table ITVl The reader should note that the num- 
bers in column 3 are not exactly the sum of the numbers 
for v — in Table UJ We have improved the accuracy of 
these by extrapolating data for all v to v — using g\{v). 

Correcting for mixing down does not change the tree 
level matching coefficients. The subtracted one loop fac- 
tors are related to the original numbers by 



a (1) -z (1) 



10 ' 



(44) 



As this subtraction is less likely to be needed in a lattice 
evaluation of NRQCD matrix elements, it has not been 
applied to the results in Tables [Til and IIIII 



TABLE IV: The mixing down subtraction. All diagrams are 
evaluated at v = 0. See the comment below Eqn. (|43[1 for 
details of column 3. 



Ma 


/; 


-^bubble H - -^earlobeH - 
-^tadpole 


-^vertex 


z 10 


4.0 


2 


0.00146 (2) 


-0.11889 (4) 


0.11743 (5) 


2.8 


2 


0.01094 (4) 


-0.17265 (6) 


0.16171 (8) 


1.95 


2 


0.03907 (5) 


-0.26196 (9) 


0.22289 (11) 


1.0 


4 


0.1870 (3) 


-0.82970 (26) 


0.6427 (4) 



B. Matrix element ratios 

If we are only interested in the ratio of leptonic widths 
of, say, T(2s) and T(ls), we do not care about the overall 
normalisation of the matrix element (which is indepen- 
dent of the mass of the meson). We can therefore express 
the ratio of leptonic widths as a ratio of differently nor- 
malised matrix elements 



^ = (Jo) + * < 
a a 



(Jo 



hi (Ji 



a 
b 2 ( J 2 ) ■ 



(45) 



The advantage of this is that for T states v 2 ~ a s ~ 0.1. 
We can obtain a ratio that is accurate to a few per cent, 
0(1% — 5%) (to two loops, effectively) by knowing a , ai 
to one loop and ai to tree level. That is, by knowing no 
more than we have already calculated in this paper: 



(o) 

= ^ _j\ a s 

U1 — 'n\ (0) 



«1 
O 



6 2 = 



02 
« 



(0) 
"0 

,(0) 



(0) (1) 



(46) 



We give these values for the unsubtracted currents in 
Tables |ITJ and IIIII Note that the inclusion of Ji at this 

as there is no mixing down 



order does not affect a 
at tree level. 



(i) 

) 



VI. SUMMARY AND CONCLUSIONS 

In this paper we have presented a method to determine 
the QCD/NRQCD matching coefficients for electromag- 
netic decays of heavy quarkonia in lattice perturbation 
theory to order 0(v 4 , a s v 2 ). This calculation was car- 
ried out for a realistic lattice NRQCD action using largely 
automated methods for performing lattice perturbation 
theory. 

The lattice NRQCD currents are given in Eq. (|34D . 
When calculating their matrix elements in a lattice 
Monte Carlo simulation, we have a choice as to whether 
we replace M by the renormalised heavy quark mass or 
the bare mass. If we choose M to be the renormalised 
heavy quark mass, the relevant matching coefficients are 
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given in Table HU If the bare mass is used instead, the 
matching coefficients include Zm and are given in Ta- 
blelnll 



We note that for the smaller quark masses, the and 

6j coefficients of the current J\ are very much smaller 
when the renormalised quark mass is used. This is partic- 
ularly relevant to NRQCD simulations of charm quarks 
on fine lattices, and shows that the use of the renor- 
malised rather than bare mass is a major source of im- 
provement in such simulations. 

Individual Feynman diagrams vary monotonically with 
the mass, but when combined together the competing de- 
pendencies lead to the final answer varying as a compli- 
cated function of M. 

We have performed a wide variety of checks of our 
calculation: we have confirmed that the Feynman rules 
are correctly generated by comparing with separately- 
obtained expressions in the literature, and that the one- 
loop self energy renormalisation similarly agrees. We 
have checked that the infrared divergences vary as ex- 
pected with changes in the size of the regulating gluon 
mass and the choice of gauge. We have also checked 
that the final answer is independent of both of these fac- 
tors. We have assured ourselves that the statistical er- 
rors quoted by VEGAS are consistent with the size of 
variations in the Monte Carlo estimates of the one-loop 
integrals. 

These results could conceivably be checked using a se- 
ries of high-/3 Monte Carlo simulations [37j • Looking fur- 
ther, the computation of the perturbative one-loop cor- 
rection to the coefficient c\ of the <r ■ B could be carried 
out using the methods employed in this paper. 



0.4 ■- 



0.2 ■ 



Q 

O 
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v aM = 4.0 
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FIG. 4: The fits to the tadpole improved data versus velocity 
dependence of (0 | Ji QQ\^°\ 

Acknowledgments 



Performance Computing Facility for advice and use of 
the SunFire supercomputer. A.H. thanks the U.K. Royal 
Society for financial support. The University of Ed- 
inburgh is supported in part by the Scottish Universi- 
ties Physics Alliance (SUPA). G.M.v.H. is supported in 
part by the Canadian Natural Sciences and Engineering 
Research Council (NSERC) and by the Government of 
Saskatchewan. 



APPENDIX A: FEYNMAN RULES 



We use an automated method to obtain Feynman rules 
from the actions and currents used in this calculation. 
The algorithm and its implementation are described in 
Ref. [30(. This allows us to specify the action as a set of 
Wilson line contours that are then Taylor expanded. The 
symmetries of the action are exploited to produce very 
compact descriptions of the reduced vertex functions as 
sums of n monomials (each involving a relatively expen- 
sive exponentiation) . 

The gluonic action expansion has been tested in a num- 
ber of calculations [U 0, H^, [4(| • The expansion of the 
currents was checked by hand. 

We tested the NRQCD action expansion by compar- 
ing with the Feynman rules quoted in Eqs. (A11-A36) 
of Ref. [l3| . We find complete agreement for general Cj , 
save in Eq. (A33) which gives the two gluon vertex for 
momenta specific to the gluon tadpole graph. Our auto- 
mated method shows this expression to be incomplete; it 
should read: 



-0.02 ■ 



-0.025 ■ 



o -0.03 ■ 



-° -0.035 - 



-0.04 



Symanzik glue 



O a=0.5 
■ a=1.0 
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Wilson glue 



-0.045 - XX 

nil i iii 

0.0001 0.001 0.01 

gluon mass, (an) 2 
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FIG. 5: Jqcd — ^nrqcd for the NRQCD action described in 
Appendix [CJ at aM = 2.1 and v = 0.03. 
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[Oi&"(k,k,q,-q) = 



f(5 



2(aM ) 3 ' An(aMo) 2 ' " 

io 2 



3 k 1 

§ij cos(^) ^ sin 2 (y ) + - sin(fci + y) sin(% + 



z=i 



16(aM ) 



-C5 



($n,j8u,o +S»,oS 1/ j) cos(fcj + y) sin(g.,) cos(y)j7 i0 



2 cos(fcj + — ) sin(9o) cos(y ) r? jQ 



12(oM ) 



cos(kj) — cos(2fcj) cos 2 (y) + - sin (2fcj) sin (qj) 



(Al) 



where the change is the addition of the final, underlined 
term. This vanishes for k = and so does not affect the 
results in Ref. 13]. Nonetheless, our detecting it high- 
lights the usefulness of an automatic action expansion 
program both for developing new improved actions and 
for checking existing perturbative results. We are happy 
to share copies of the program with interested parties. 

The NRQCD action in Eq. naturally factorises 
into the product of several distinct operators: 

5 nrqcd =^1^-4 AtBtUlBt-lA^ , (A2) 



where 



A^l- 
B=U- 



2n J 
aSH 



(A3) 



and the subscript refers to the timeslice on which the 
fields are located. 

In the "by-hand" expansion it simplifies the algebra 
to derive separate Feynman rules for A and B and com- 
bine them using the convolution theorem (l3l . HHj. We 
also follow this approach: the AB and BA factors are on 
different timeslices so no compression of the set of mono- 
mial factors ( "entities" ) contributing to the reduced ver- 
tex function can occur. Without such compression, it is 
computationally cheaper to calculate the Feynman rules 
as a convolution of the expansions of A and B. The im- 
plementation of this has been checked by comparing with 
the reduced vertex functions from the expansion of the 
full action. 

Partial derivatives of the Feynman rules are computed 
automatically in the code as per Ref. [30( . We exploit the 
fact that the velocity is purely along the z-axis to write 

" 1 8 (A4) 



dp 2 2p 3 dp 3 
The total on-shcll derivative is implemented as 
d d dpo d 

dp 2 dp 2 dp 2 dpo 
dpo i dT{p) 

dp 2 l-T(p) dp 2 ' 



with T(p) = G o 1 (0,p) coming from the bare fermion 
propagator. 

APPENDIX B: RENORMALISING THE 
FERMION PROPAGATOR 

In this appendix we review the one loop renormalisa- 
tion of the fermion propagator. The bare fermion prop- 
agator is 

aG- 1 (p Q ,p) = l-z(l- aT{p 2 )) (Bl) 

where z — e~ zapo and T(p) is the kinetic energy. The 
0(a s ) NRQCD quark self-energy can always be written 
as 

a£(p ,p) = A + B{ Po ,p)aT(p) + 

C(p ,p)[l-z(l-aT(p))} , (B2) 

where A is a constant. The resummed propagator is 

aG~ 1 (p ,p) = aGQ 1 (p ,p) - a s aY,(p ,p) 

= (1 - a s (A + C)) [1-2(1 + a s A) x 
(1 - aT{p) [1 - a s B/z])] 
+0{a 2 s ) (B3) 

In the infrared limit of small p 2 , this should be compared 
to the renormalised form of Eq. (|B1[) : 



aG-^Z-^l-zll-aTnip)}) , (B4) 

with T R (p) = p 2 /{Z M M). Identifying z = z(l + a s A), 
the additive shift in the rest energy is 

aAE rcst = ln(z/z) 

= -a s A + 0{a 2 s ) , (B5) 

and A — aT,(p a — 0,p = 0). The po pole in the propa- 
gator occurs at z — z = (1 — ciTr) . The wavefunction 
renormalisation is found by Taylor expanding Eq. (|B3D 
around this pole: 



dB 

ZAp) = l + a s [A + C + BaT-aTz— 

oz 



(A5) 



= 1 + a s \ aS + 



d{ia Po ) J on . shell 



on-shcll 

(B6) 
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where the expressions are evaluated on the mass shell. As 
the terms in brackets are already 0(a s ), it is sufficient to 
identify T and Tr and evaluate them at the pole of the 
bare propagator z = (1 — aT)^ 1 . This result is general 
and includes all orders in v 2 at 0(g 2 ). Morningstar [28| 
gives the expression for Z$ at zeroth order in v 2 and our 
result agrees with his to this order. 

Working on the renormalised mass shell, the mass 
renormalisation follows from 



1 



dT R 



2Z M M dp 2 



p 2 =0 



1 

2M 



d(BT) 



dp 2 



(B7) 



p 2 =0 



We note that the total differential must also be evaluated 
on the (bare) mass shell. From this we obtain 



1 + a s 2M 



daT, 



dp 2 



(B8) 



Tadpole improvement affects Ai5 ros t and Zm, but not 
Z^p for NRQCD actions that arc symmetric under time 
reversal Q. The one-loop contributions are given in 
Eqs. (35, 36) of Ref. 



APPENDIX C: FURTHER CODE TESTS 

In this appendix we describe further, non-trivial tests 
of our perturbative calculation that verify that our con- 
tour shifting and numerical integration techniques are 
correct. 



1. One-loop self energy 

We have calculated the renormalisation of the NRQCD 
propagator as per Appendix [B] using Feynman gauge and 
a gluon mass of {an) 2 = 1CT 4 . The results are given m 
Table [Vj For comparison, we also give the results of 
Gulez et al. [HI]. Our data agree very closely. This 
provides further evidence that not only are our Feyn- 
man rules correct, but also that we are combining them 
correctly to form diagrams and evaluating the resulting 
integrals correctly using VEGAS. We have also checked 
that the results are correctly gauge variant and that the 
effect of the finite gluon mass is negligible. 



2. Gauge covariance and invariance 

We have also looked closely at the effect of changing 
the gauge and infrared regulator. For these tests, we 
use a simpler NRQCD action with coefficients c* = 
for i — 1 . . .4 and C5 = cq = 1, as used in Ref. [r|, 
with n = 2 and aM = 2.1. We used both the Wilson 
and Symanzik-improved gauge actions and set current 



J 1 = 0, which implies I c . 



rlobc 



'bubble 



'tadpole 



0. 



We use three choices of gauge: Feynman gauge A = 1, 



an unnamed gauge with A — 2 and Yennie-Fried gauge 

Firstly, Iqcd — ^nrqcd should be independent both of 
the choice of gauge and the gluon mass. This is seen for 
v = 0.03 in Fig.[5l We note that the size of the scatter of 
points about a single mean value is consistent with the 
statistical errors assigned to the data points by VEGAS. 
This gives us some confidence that these errors are not 
being underestimated in our calculation. 

Next, Z and V separately have infrared divergences 
that are regulated by the gluon mass, but which cancel 
in Z+V. The cancellation has already been shown by 
the absence of diverging behaviour at small a/i in Fig. [5] 

Here we check that the individual diagrams show the 
correct divergence. We compare with the expectations 
for continuum QCD: lattice NRQCD is an effective de- 
scription of this and must preserve the same infrared 
structure (up to possible discretisation errors of order 
an). 

Given the lack of overall divergence in Z + V, it is 
sufficient to concentrate on Z, which is determined to 
greater statistical accuracy. 

The one- loop continuum expression is given in Eq. (7- 
44) of Ref. [If (adding a colour factor of 4/3): 



Z r 



1 



12tt 2 



1 A 2 

A M 2 



ft 



The infrared divergent contribution is 



9 
4 
(CI) 



Zm — — - 



1 



12tt 2 



3- 



ln/r 2 



(C2) 



which vanishes in Yennie-Fried gauge. 

In Fig. [S] we plot Z with the continuum divergence 
removed (replacing (i by a/i). There is no discernible 



-0.02,£K> 



-0.04 ■ 



o 



-0.06 ■ 



-0.08„ L 



O a=0.5 
■ a=1 

A a=3 



0.001 0.002 
gluon mass, (a\i) 2 



0.003 



FIG. 6: Z - Z m , for the NRQCD action described in Ap- 
pendix [C"2l at aM = 2.1 and v = 0.03 using the Wilson gauge 
action. 
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TABLE V: The renormalisation of the fermion propagator, as compared to Ref. [l3l ]. Note that an IR factor corresponding to 
Eq. (|C2[I with (afi) 2 — 1CP 4 has been applied to the subtracted data quoted in Ref. [l^] for Z^\p = 0). 

4 1} (P = 0) Z§ aAE™ 



aM 


n 


Us 


Ref. [13] 


Us 


Ref. [13] 


Us 


Ref. [13] 


4.0 


2 


1.8207 (7) 


1.813 (3) 


0.0817 (5) 


0.082 (4) 


0.8390 (1) 


0.850 


2.8 


2 


1.6232 (6) 


1.617 (3) 


0.2350 (6) 


0.235 (4) 


0.7570 (10) 


0.767 


1.95 


2 


1.3494 (5) 


1.344 (3) 


0.4201 (8) 


0.421 (4) 


0.6765 (10) 


0.689 


1.0 


4 


0.4334 (7) 




0.8285 (16) 




0.9684 (13) 






6 




0.410 (3) 




0.859 (4) 




0.758 



divergence as a/i — > 0. The slight gradient betrays a 
residual dependence on the gluon mass. To emphasize 
this, we plot the deviation A2(o/j) = Z(a/j.) — Z(10a[i) 
in Fig. [71 The deviation disappears as we take a/i to zero 
and is a discretisation effect. 



3. Current matching at v = 

Finally, we have tried to verify the one- loop, O(v ) 
annihilation current matching of Jones and Woloshyn 
[Tfj| . Following the method in the main text, we get 
a W = _o.0225 (3) for (a/i) 2 = 10~ 4 and -0.0228 (3) 
for (a/i) 2 = 10~ 3 (using n = 2 and aM = 2.1 with the 
Symanzik gauge action). The extrapolations to v = 
are shown in Fig. [5] The statistical compatibility of the 
results shows that (a/i) 2 = 10~ 4 is small enough that any 
residual gluon mass dependence of the results is swamped 
by the statistical uncertainties in the VEGAS integration. 

At the same parameter values, Jones and Woloshyn 
give a,Q = —0.0253 (3) [inserting the appropriate num- 
ber from Table II into their Eq. (27)]. This result broadly 
agrees with ours, which gives us confidence that there are 
no gross disagreements in our methods. 



There is still a small, but apparently significant devi- 
ation, which we also see at a second mass value. The 
stringent tests described in these Appendices were our 
attempt to account for this difference. As already de- 
scribed, we have checked our Feynman rules are correct 
and give the correct self energy (and derivatives). We 
find the correct infrared divergences and Lorentz invari- 
ance. We have gauge invariance and independence on 
the gluon mass regulator. We have also checked that the 
statistical errors quoted by VEGAS are not underesti- 
mated. In the light of these, we feel confident that our 
calculation is correct. 
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FIG. 7: AZ for the NRQCD action described in AppendixlC2l 
at aM = 2.1 and v = 0.03 using the Wilson gauge action. 
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FIG. 8: Determination of aj, for the NRQCD action de- 
scribed in Appendix IC 2l at aM — 2.1 and v = 0.03 using the 
Wilson gauge action. 
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